Quantum annealing of an Ising spin-glass by Green's function Monte Carlo 
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We present an implementation of Quantum Annealing (QA) via lattice Green's function Monte 
Carlo (GFMC), focusing on its application to the Ising spin-glass in transverse field. In particular, 
we study whether or not such method is more effective than the Path-Integral Monte Carlo (PIMC) 
based QA, as well as classical simulated annealing (CA), previously tested on the same optimization 
problem. We identify the issue of importance sampling, i.e., the necessity of possessing reasonably 
good (variational) trial wavefunctions, as the key point of the algorithm. We have considered two 
possible classes of trial wavefunctions, a mean-field single-site one — whose optimization is however 
a very difficult task — and a Boltzmann-like choice. We performed GFMC-QA runs using such 
a Boltzmann-like trial wavefunction, finding results for the residual energies that are qualitatively 
similar to those of CA (but at a much larger computational cost), and definitely worse than PIMC- 
QA. We conclude that, at present, without a serious effort in constructing reliable importance 
sampling variational wavefunctions for a quantum glass, GFMC-QA is not a true competitor of 
PIMC-QA. 

PACS numbers: 03.67.Lx, 75.10.Nr, 03.65.Xp, 02.70.Ss, 05.10.Ln, 07.05.Tp 



I. INTRODUCTION 

Quantum annealing (QA) 1 is based on the idea of 
searching for the ground state of some classical Hamilto- 
nian by adiabatically switching-off an appropriate source 
of quantum fluctuations, in much the same way as tem- 
perature would do in thermal annealing. This approach 
is also known as adiabatic Quantum Computation, in the 
Quantum Computing community^. For recent reviews of 
the QA field, see Refs. Qj. 

A very popular QA approach is based on an imaginary- 
time Quantum Monte Carlo (QMC) implementation, i.e., 
the Path-Integral Monte Carlo (PIMC) approach. A cer- 
tain success in the application of PIMC-QA has been 
obtained in most of the cases studied: The folding of 
off-lattice polymer models^^ the random Ising models 
and the random-field Ising model ground state search^, 
Lennard-Jones clusters optimization^iS and the travel- 
ing salesman problem— Nevertheless, a counterexam- 
ple exists^ where PIMC-QA performs definitely worse 
than simple CA: The 3-Boolean-Satisfiability (3-SAT) 
problemp 13 which is a prototype of a large class of 
hard combinatorial optimization problem (the so-called 
nondeterministic polynomial- (NP)- complete class, see 
Ref.0). 

Generally speaking, the PIMC-QA failure depends on 
the properties of the "landscape" of the problem at hand 
or even to a bad performance of the implementation. 
In order to understand its features in details, in a re- 
cent paper— we have studied the PIMC-QA performance 
focusing our attention on a simple, but highly instruc- 
tive, toy-problem: The double-well potential. There we 
learned a few potential dangers of the PIMC-QA method, 
in particular: i) The unavoidably finite temperature T 
used in the simulation, which provides a thermal lower 



limit to the average residual energies attained by the 
algorithm, ii) The possible severe difficulties (ergodic- 
ity breaking) in sampling the PIMC action close to a 
Landau-Zener crossing of ground state levels. 

We propose here to investigate a different Quantum 
Monte Carlo (QMC) based QA algorithm, as an alter- 
native to PIMC-QA. A very natural choice is provided 
by the Green's Function Monte Carlo (GFMC). GFMC 
is different from PIMC since it can directly sample the 
ground-state of a quantum Hamiltonian, avoiding, in 
principle, one of the PIMC drawback, i.e., the finite tem- 
perature T. However, contrary to PIMC, GFMC de- 
mands the knowledge of good variational wavefunctions 
to implement what is called the importance sampling. We 
will appreciate soon how serious a drawback can this be. 

A natural benchmark for a test of this new GFMC-QA 
algorithm is provided by the random Ising model, a chal- 
lenging optimization problem already addressed through 
PIMC-QA, as well as standard classical simulated an- 
nealing (CA), in the recent past&i The Hamiltonian of 
random Ising model in transverse field is: 



H(T) = -£ Jij ola] -Tj2°i=Hc 



kin 



(i) 



where j) indicates a sum over nearest-neighbors, J^j 
are random nearest-neighbor Ising coupling constants, 
and of , of are Pauli's matrices on lattice site i. If we de- 
note by {Si} a generic configuration in the Hilbert space 
(where Si = ±1 are the eigenvalues of of matrix), the 
classical function we want to minimize is just given by 
the first term in Eq.|H i.e., E cl ({S t }) = ({Si}\H cl \{Si}) 
- the Ising glass energy itself — , which plays the role 
of potential energy. The second term in Eq. ^ -Hfcm = 
— r of, is the needed source of quantum fluctuation, 
which plays therefore the role of a kinetic energy. More 
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in detail, we will concentrate our efforts on a problem 
instance which has been analyzed extensively in Ref. 00: 
it is a two-dimensional (2D) instance, on an L x L square 
lattice with L — 80, and the couplings Jij drawn from a 
flat distribution in [—2,2]. 

Since we want to employ Eq. ^ as the Hamiltonian of 
our QA dynamics, the transverse field T represents the 
annealing parameter of the system; the goal is to follow 
the time-dependent dynamics with a F(t) which starts 
from very large values, and is ramped down to zero in 
a certain annealing time r. We emphasize that such a 
transverse field term is not just a theoretical concept. 
Indeed, the whole field of QA was strongly revived by 
experimental results on the disordered Ising ferromagnet 
LiHoo.44Yo.56F4, where the transverse external magnetic 
field r was actually applied to the system and manipu- 
lated in the laboratory, to perform a true Quantum An- 
nealing experiments^* 

The rest of the paper is organized as follows. In Scc.lTTl 
we present the main ideas of a GFMC-based QA ap- 
proach, with a sketch of the main ingredients of the al- 
gorithm, for the benefit of the non-expert reader. In 
Sec lIIII we present the results of our variational study 
of trial wavefunctions, showing the inherent difficulties 
associated to the selection of good wavefunctions in a 
disordered quantum system. In Sec II VI we present the re- 
sults of fixed-r GFMC calculations. In Sec[V]we discuss 
the GFMC-QA results, and compare them with previous 
PIMC-QA and CA data on the same problem. Finally, 
in Sec. I VII we give some concluding remarks. 

II. GREEN'S FUNCTION MONTE CARLO 
QUANTUM ANNEALING: IDEAS 

The ideal scope of a QA approach is to take some initial 
state V'(O) and let it evolve, according to the Schrodinger 
dynamics associated to a time-dependent Hamiltonian 
H(t) interpolating between an extreme quantum regime 
and the classical problem one is interested in: 

i^ = H(t)4>(t)> (2) 

where we have set h = 1. For instance, for the prob- 
lem we have set to consider, we could take H{t) — 
- £«j) •/..., ofoj ~ r(t) J2 t of, where T(t) is initially 
very large, and slowly decreased towards zero, as illus- 
trated in Fig. ^ As argued in Ref. [IH, an imaginary- 
time Schrodinger evolution would be, for optimization 
purposes, equally good, and most likely even superior to 
the standard real-time evolution. With this in mind, we 
can get rid of the i in the time-derivative term in Eq. [21 
substituting it with the imaginary time prescription — dt- 
If we also imagine the gradual decrease of T to be made 
step- wise, as sketched in Fig.^ then the solution for tp(t) 
is obtained by repeated applications of an imaginary-time 
propagator 

V>(r) =e- H ^^---e- H ^ T ^(0) , (3) 



where r = n + ■ • • r n is the total annealing time, Ti > 
r 2 > • • • > T„ ~ is a decreasing sequence of transverse 
fields, and H(Ti) is a shorthand for H(t) with a value Ti 
of the transverse field. Each application of the imaginary- 

r(t)' 1 




FIG. 1: Sketch of the annealing of the transverse field F, 
linearly in a time r = Ti+T2 + - • • T n , or in a step- wise fashion. 

time propagator e~ H ( ri ) Ti effectively tends to filter out 
the corresponding ground state of H(Ti) from the state to 
which it is applied. In turn, e~ Ht ^ r ^ Ti can be obtained as 
a repeated application of many infinitesimal propagators 
of the type [1 - AtH(Ti)], i.e., 

V(t + n) = e- H ^m = H [1 - AtHpi)] m ■ (4) 

The Green's function Monte Carlo (GFMC) is just a 
stochastic technique which implements such a form of 
propagation. More precisely, if we define, recursively, 

^ n+1 (x') = Y / G i pJ Jn (x) , (5) 

X 

where ip n (x) = (x\ip n ), \x) being a shorthand for a generic 
spin-configuration describing the Hilbert space of the 
problem, and 

Gi% = (x'\G^\x) = ( x i\l~AtlH(T)-E T }\x) 

= (l + AtE T )S x ,, x -At(x'\H{T)\x), (6) 

one can show that - for large n - the iterated state tp n 
converges (apart from a normalization constant) to the 

ground-state iPqs( x ) °^ ^ ^ ^ s chosen to be suit- 

ably smalls (Et is an estimate of the ground-state en- 
ergy which allows to reduce the statistical fluctuations, 
see also Ref.l2fjh 

The problem in Eq. [S] looks superficially similar to an 
ordinary Markovian Master equation, with a few very 
crucial differences: i) The ip n {x) are not probabilities, 

but amplitudes] ii) The Green's function in Eq. [SI 

unlike the transition probability of a Master equation, 
is not necessarily made of non-negative elements, and 



3 



is, in general, not column- normalized, ~52 x iG x , x ^ ^ 

unlike a Markov transition probability. In summary, the 

process underlying the iterated-power method is not a 

properly defined Markov chain, and, therefore, it cannot 

be immediately simulated, as it stands, with a Monte 

Carlo approach. 

Problem ii) above can be quite serious: If some of the 
(r) 

matrix elements of G x , x are negative, no possible inter- 
pretation of it as a "transition probability" is possible. 
This is at the heart of the so-called sign problem*^ in 
Quantum Monte Carlo. In the following, we will assume 
that a choice of basis is possible in which no sign prob- 
lem exist, i.e. all matrix elements of G are non-negative, 
G x ',x > 0. This is certainly true for the Ising glass in a 
transverse field. More generally, since the choice of the 
kinetic energy to be used in QA is at our disposal, it is 
wise to choose the signs of Hkin in Eq. ^ such that no 
sign problem occurs. Still, we miss the correct column- 
normalization, 

E G % = l + &tE T -AtJ2 H x -, x = f b x + l, (7) 

x' x' 

A way out of this difficulty is to factorize G^ T ' in terms 
of a stochastic matrix p x > , x ~ by definition, a matrix with 
all positive elements p x ',x > 0, and with the normaliza- 
tion condition ^2 X , p x ',x — 1 for all columns - times the 
scale factor b x defined above. Indeed, with the previous 
definition Q of b x , the matrix 

= G x T ,]jb x . (8) 

is trivially positive and column normalized and, there- 
fore, it is a suitable transition matrix for a Markov chain 
in x-space. 

The crucial idea is then to extend the configuration 
space where the Markov process is defined, adding to the 
x a non-negative weight factor (hereafter, the weight) 
w. This extended configuration space will be labeled by 
(x,w). The pair (x,w) is often called a walker, because 
it will be the basic entity in the Markov chain "random 
walk". The weight part will take care of b x , while x 
will be taken care of by p x > ,x- More precisely, if (x„, w n ) 
indicates a walker at iteration time n, in this extended 
configuration space, we set up the following Markov pro- 
cess: 

a) generate 

^n+i — x 1 with probability p x ',x n 

b) update the weight with w n +i = w„b x ■ (9) 

In words: The walker performs a random walk in the 
Hilbert space x of the system and in the weight space w; 
such a random walk is composed of a standard Markov 
chain in x-space, associated to the p x ',x, plus a multi- 
plicative process for the weight w n — > w n +i = w n b x . 
By moving in this way, the walkers visits every point 
in the (x,w)-space with a probability P n (x n ,w n ) whose 
first moment can be shown to be proportional to the 

lf> n (x) OC / dw n W n P n (x,W n )m± 



Eq. [5] is the basic version of a GFMC algorithm^ 
In this form, however, the algorithm simply does not 
work in practice. The reason for this failure is not dif- 
ficult to grasp. While x n — ► x n +\ is a plain Markov 
process, the weight update w n — > w n +\ — w n b Xn is a 
multiplicative process with random factors b Xn , which 
is prone to very large fluctuations>i2iSi w n might grow 
large, or become negligibly small, in just a few itera- 
tions, and the whole algorithm would go wild, because 
error bars in the calculations of the averages grow in an 
uncontrolled way. The cure to this disease goes through 
the introduction of many walkers and through perform- 
ing occasional "reconfigurations" of their weights, via the 
so-called branching^ In practice, one propagates simul- 
taneously a set of M walkers defined by weights Wi and 
configurations Xj, for i — 1, • ■ ■ M. Before the variance of 
the weights Wi becomes too large, one appropriately re- 
defines the set of walkers - by reproducing some of them 
and deleting some others - in such a way as to drop those 
with excessively small weight, and to generate copies of 
the more important ones^ 

The analogy of such a many-walker GFMC with a 
genetic-like algorithm is worth noting. Each walker 
(x, w) plays the role of an individual that propagates 
(mutates) increasing or decreasing its fitness - repre- 
sented by the accumulated weight w, related to the wave- 
function amplitude ip( x )- A mutation is here simply a 
step of the algorithm, which attempts a single spin-flip 
of a certain site in the configuration: this is what the off- 
diagonal matrix elements (xi\H\x) do. At certain times, 
branching occurs, which modifies the population of in- 
dividuals by favoring the survival of those with highest 
fitness (largest w). The only genetic feature that is miss- 
ing in the quantum mechanical case is the possibility of 
cross-breeding (mixing of genetic codes of two configu- 
rations, to give rise to new configurations): this would 
correspond to non-local moves which change the config- 
urations in a global way. 

The final, important, ingredient that makes the algo- 
rithm work is the so-called importance sampling^ It can 
be seen, in the genetic analogy proposed before, as a 
way of proposing mutations (single spin flips) that in- 
stead of being equally probable, with a matrix element 
F, are biased by a function which guides the system to- 
wards the most representative configurations. More pre- 
cisely, suppose we have a reasonable guess of the ground 
state i^gs (x) in the form of some nodeless wavefunction 
^t(x), known as trial (or guiding) wavefunction. It is 
then enough to substitute G with the so-called impor- 
tance sampling Green's function: 

G x ,, x =iPt(x')G x ,, x ^ 1 (x) , (10) 

which just rescales by an extra factor iPt{x')/i^t{x) a 
transition from x to x , thus favoring those transitions 
where iPt(x')/iPt(x) is largest. In general, G x ^ x is not 
symmetric, but one can still apply to it the same decom- 
position in defining the corresponding Markov chain 
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© with: 

Px' ,x — G x ' ,sc / &x 



K. = 



\G X , X = l + AtE 7 



At 



= I -At (E L (x)-E T ) , (11) 
where the local energy El(x) is defined as: 

(ih\H\x) 



E L {x) 



{iPt\x) 



(12) 



If the guessed trial wavefunction -0T coincides with the 
actual ground-state wavefunction iPt(x) — iJjgs(x), then 
El{x) = Eqs is a constant, and one can show that statis- 
tical fluctuations in the calculation vanish exactly. This 
is the so-called zero variance property^ Therefore, by 
variationally improving the quality of the guiding func- 
tion iI>t{x) one can substantially reduce the error bars in 
the calculation. 
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FIG. 2: (Color online) The average best residual energy ob- 
tained by GFMC-QA for the 80 x 80 instance of the random 
Ising model studied in RefsUl versus the total annealing- 
time t. Upper rhombi: GFMC-QA results without impor- 
tance sampling (ipT = 1). Lower rhombi: GFMC-QA results 
with importance sampling performed by using the optimal 
trial wavefunction if>!P of Sec. Ell The GFMC time-unit 
is a single spin-flip, while CA and PIMC-QA Monte Carlo 
time units are sweeps of the entire lattice (see Ref. 0). The 
transverse field is linearly reduced down to 10~ 4 in a total 
annealing-time r, starting from Fo = 2.5. We used here 
M = 20 walkers and performed branching at every MCS 
(ts = 1). Previous results obtained by classical simulated 
annealing (CA) and by Path-Integral Monte Carlo quantum 
annealing (PIMC-QA) with P = 20 Trotter slices^ are shown 
for comparison. 

The fact that importance sampling is indeed a crucial 
ingredient is demonstrated, for our case, by Fig.|2 where 
we show the results obtained by GFMC annealing with- 
out importance sampling (top curve) compared with the 



results obtained with importance sampling, which we will 
illustrate later on. Quite evidently, the residual energy 
obtained without importance sampling is terribly bad: 
any short classical simulated annealing would do better. 

Finding a reasonable trial function iPt(x) for the prob- 
lem at hand is therefore an essential part of a GFMC- 
QA application, and constitutes the delicate point of the 
whole algorithm. In the next section we will describe the 
choices of ipr(x) we have tested for the Ising case, and 
the difficulties encountered. 



III. VARIATIONAL WAVEFUNCTIONS FOR 
THE ISING SPIN-GLASS 

Finding a good trial wavefunction for a random Ising 
model in a transverse field is a highly non-trivial task. 
The first idea that comes to mind is a kind of "mean- 
field" wavefunction, made up of a product of single-site 
factors as: 



\4 MF) ) 



N 

n 

i=l 



i)i 



y/2 cosh(hi) 



(13) 



where {hi}, the local fields on each site i, are varia- 
tional parameters to be optimized for each given value 
of the transverse field T. The optimization of the {hi} 
amounts to finding the minimum of the variational en- 

i.j) m i m i 



ergy, E y T 



where = tanh(ft,i) are the local magnetizations. The 
stationarity conditions required by the minimization, 
d E^ IF ^ jd hi = 0, read for each site i: 




Ji j rrij I + r mi 



(14) 

where N(i) indicates the set of nearest-neighbors of site 
i. As it turns out, finding solutions of Eq. (|14fl with opti- 
mal variational energies is simple only for large enough T, 
where the quantum paramagnetic solution hi = mi = 
is found. Such a solution, representing the T = +oo 
ground state with all spins aligned along the +x direc- 
tion, survives down to some value T cr of the transverse 
field, below which non-trivial solutions of Ea. 1141 - with 
non- vanishing local magnetizations rrij ^ - start to ap- 
pear. However, the number of those solutions (local min- 
ima) is large. In the low-r region, our minimization prob- 
lem is just the quantum counterpart of the well-known 
Weiss mean-field approach for the classical random Ising 
model, 24 which is known to run into difficulties in the 
classical glassy phase. In a sense, minimizing E^ 1F ^ in 
the low-r glassy phase is not much simpler than find- 
ing the classical ground state of the problem for T = 0: 
We transformed a minimization task in a discrete space 
of variables, Si = ±1, into one where the variables are 
continuous, to^ G (—1, 1), but the task itself is of com- 
parable difficulty. We anticipate that difficulties similar 
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FIG. 3: (Color online) Distribution of the variational energy 

per spin eimr = nii n {fe-}{V'T ^'I-HI^t F ^}/N obtained by 
optimizing, with Conjugate Gradients, the local-fields hi in 
the trial wavefunction \ip'^ F ^) defined in Eg. 1131 These re- 
sults refer to the case of a fixed F = 0.1. The upper panel 
has been obtained by repeated conjugate gradient (CG) mini- 
mizations, while in the lower panel we have used an annealing 
CG scheme (see text). In both cases, the initial local fields hi 
are randomly taken from the interval (—0.5, +0.5). 



to the classical Weiss approach will plague our search for 
the minima of 12^ in the low-r phase. 

Given the explicit analytic expression for the energy 
to be minimized, we made use of a standard Con- 
jugate Gradients (CG) algorithm to find local min- 
ima solutions for the {hi}, and the corresponding op- 
timal value of the variational energy per spin, e^ r F ^ = 
mm{ h .}ET IF \{hi})/N. However, as anticipated, the re- 
sults obtained for ei^- when T < T cr depend on the ini- 
tial conditions of the algorithm, so that the only meaning- 
ful thing to do is to show histograms of e^ F ^ , obtained by 
repeated - and uncorrelated - CG- minimization searches. 
In Fig. |3 we report the results of two different calcula- 
tions for r = 0.1, which turns out to be in the low-r 
glassy phase. The top panel displays the histogram of 
500 repeated searches by CG-minimization, each search 
starting from initial hi which are randomly distributed 
in the interval (—0.5, +0.5). The large solid arrow marks 
the location of the classical ground-state energy per spin, 
£gs = Eqs/N ~ —1.5805167, as a reference. In most of 



the attempts, we simply find minima with e^ F%> ~ 0, 
that have nothing to do with the classical GS which 
should dominate for such low T value; only rarely, we 
end up with local minima in the range ei^r — —1.215, 
which is still quite far from ecs- Although, in princi- 
ple, quantum effects due to the finite value of T make 
e^ IF ^ ^ e GS , the small T value used (T — 0.1) does 
not justify such large difference between the two quan- 
tities. This is a clear suggestion that the bare CG- 
minimization approach applied to the single-site varia- 
tional Ansatz £^ MF ' is prone to get stuck in the high- 
energy paramagnetic phase. Data (not shown, but mon- 
itored during the simulation) on the average magnetiza- 
tion, the Edward- Anderson parameter^ 5 and the super- 
position integral between the optimized variational state 
and the true ground-state, are also in agreement with 
this hypothesis. To give further support to this conclu- 
sion, we show in the bottom panel of Fig. [3] the results 
obtained by an improvement of the CG scheme, which 
is a sort of "variational" annealing. There we started 
from {hi} random distributed in (—0.5,0.5) and we per- 
formed a first CG-minimization at To = 4, a sufficiently 
high value which provides an initial paramagnetic state. 
Then, the transverse field T was repeatedly decreased 
in 10 steps of AT = 0.39 down to the desired value of 
r = 0.1, and at each step, the new minimum was searched 
via CG-minimization starting from the optimal solution 
found in the previous step. This "variational annealing" 
scheme provides results for e>^ IF ^ which are distributed 
much closer to the classical ground-state energy, cgSj al- 
though we still obtain, in a non-negligible fraction of the 
attempts, paramagnetic solutions. 

Summarizing, using the single-site trial wavefunction, 
Eq. we find that at small r w 0.1 the system is very 
unlikely to stay in its magnetized phase (i.e., the chal- 
lenging glassy state). The lack of regularity in the vari- 
ational minimization, consequence of the multiplicity of 
metastable minima, makes t/j^^ not a particularly suit- 
able guiding wavefunction for a GFMC simulation. 

A second, quite natural, choice of trial wavefunction is 
a Boltzmann-like wavefunction of the form: 



^\{Si})=N{P)e^W E ^ s ^ 



(15) 



where 1/(3 plays the role of an effective temperature, with 
[3 a variational parameter to be optimized, and E c i({Si}) 
(see Eq. is the classical energy of a given configuration 
{Si}. N{(3) is an appropriate normalization factor, which 
we will not need to calculate. Once again, for large T we 
expect to find [3 — (the exact T = +oo solution), while, 
by decreasing T, larger and larger values of (3 will favor 
regions where the "potential energy" E c i({Si}) has a lo- 
cal minimum, until we get, for r = 0, to the asymptotic 
limit (3 — > oo (ideally), required by a wavefunction which 
is perfectly localized in the global minimum (see below 
for a discussion of this point). 

To calculate the expectation value of energy with the 
Boltzmann-like choice in Eq. 1151 as a function of the sin- 
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FIG. 4: (Color online) (Top) Plot of the optimal /3, op t, 
for the "Boltzmann" trial wavefunction |V"t ) defined in 
Eq. I15I for several values of V. The dashed line is a 
guide to the eye. (Center) Optimal variational energies 

4ot° ltz) = (^T° pt) \ H \i>T° pt) )/ N corresponding to the /3 0pt 
shown in the Top panel, and the GFMC estimate of the 
total energy per spin {H}/N, versus F. The inset mag- 
nifies the small-r region, where small differences are no- 
ticeable. (Bottom) The variation residual diagonal energy 
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l3 opt shown in the Top panel, and two GFMC estimators of the 
residual diagonal energy e res = {H c i)/N — tGs- The mixed av- 
erage, and the Ceperley correction (see text). The dashed hor- 
izontal line represents the best residual energy ever achieved, 
for r > 0.01, by employing the mean-field trial wavefunction 
in Eq.[T31 



gle parameter (3, we used a standard Variational Monte 
Carlo (VMC) algorithm. — Fig. ^ shows (top panel) the 



r^(Boltz) 

energy E\ ot 



(4>T 



(.0) 



H\i/>. 



(40 \ 



for several values of the 



transverse held T. Notice that [3 op t saturates for small V 
to about (3 pt ~ 2, somewhat surprisingly at first sight, 
since, for V — > 0, one would expect (3 op t — * +oo, in such a 
way that the classical ground state dominates (i.e. Eg. 1151 
becomes a delta-like function localized in the exact clas- 
sical ground-state). This is clearly an effect of severe er- 
godicity loss of the VMC algorithm, which is not difficult 
to understand. For a given /3, indeed, the VMC samples 
a probability distribution \ifx^' {x)\ 2 = e^ /3Ec '^ with sin- 
gle spin-flip moves: its efficiency in exploring the phase 
space, therefore, is exactly identical to that of a classical 
Metropolis Monte Carlo at the temperature T = 1//3. 
Finding an optimal (3 for a given T is therefore totally 
equivalent to ask what is the effective temperature of a 
classical Ising spin glass which provides the best approx- 
imation to the wavefunction of a quantum Ising glass at 
zero temperature and non-zero T. Now, from classical 
spin-glass physics^S& we know that a threshold energy 
E t h exists below which the system has a finite complexity, 
i.e., it displays an exponentially (~ exp V) large number 
of metastable minima. Close to this threshold energy, 
the relaxation of any local algorithm towards equilib- 
rium becomes exceedingly slow (the algorithm gets stuck 
for a long time in every minimum visited) and the av- 
erage quantities measured are not representative of their 
true thermodynamical values. Evidently, for r — > 0, the 
variational algorithm is not visiting the regions near the 
true minima of the classical energy, but is wandering in 
a high-energy band of metastables states, separated by 
moderate energy barriers. In such small and 

finite value of (3 allows to still overcome such barriers, 
so as to find slightly more favorable local minima, while 
perfect localization {[3 — > +oo) in a wrong exited state 
would lead to an average bigger residual energy, 

The central and bottom panels in Fig.0|show the opti- 



mal variational energies e^°'* z ^ 



(Boltz) 
Cres 



)/N, 



and the variational residual energy 
{^° pt) \H cl \^° pt) )/N - e GS corresponding to the 
optimal [3 shown in the top panel, for several values of 
transverse field T. For large T values, the variational 
total energy (center panel) is linear in T, as it should, 
since the transverse field kinetic term dominates in the 
quantum paramagnetic phase (see Eq. while the 
variational residual energy per-site is of order 1. By 
decreasing T, we notice that the variational residual 
energy saturates, for small T, to finite non-zero values, of 
order 0.03, in agreement with the previously noted sat- 
uration in the optimal f3 op t, due to ergodicity breaking. 
A closer inspection shows that the variational residual 
energy is actually non monotonic for T < 0.25, again an 
artifact of sampling difficulties. Notice, however, that 
this saturation value is definitely below the best (down 
to r = 0.01) results provided by the previously discussed 
^(mf) ^ wn j cn j g £ orc j er 0.035 (shown for comparison by 
a dashed horizontal line). Therefore, with all its pitfalls, 
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the Boltzmann-like trial wavefunction in Eg. 1 151 provides, 
at low r, a marginally better approximation of the true 
GS, than that obtained by the mean-field Ansatz, 
Moreover, ip^ is a l so much better behaved, and 
simpler to optimize, as far as the minimization problem 
is concerned. For these reasons, we have decided to 
work out our GFMC results using the Boltzmann-like 
wavefunction only. 

IV. BEYOND VARIATIONAL: GFMC AT 
FIXED TRANSVERSE FIELD 



the mixed average estimates, which, in turn, are lower 
than the variational result. However, the non-monotonic 
behavior of the residual energy data for small T < 0.25, 
previously noted for the variational results, should ring a 
bell about the quality of the trial wavefunction, and the 
efficiency of the sampling (i.e. ergodicity), in that region. 

Summarizing, we have shown that a fixed-r GFMC al- 
gorithm, while improving a bit on the variational results, 
still suffers from the rather unsatisfactory quality of the 
importance sampling trial wave- function, and still shows 
ergodicity problems in the small-r region. 



The variational results presented in the previous sec- 
tion can, in principle, be improved by GFMC. Before 
showing the actual GFMC annealing results, it is in- 
structive to see the performance of GFMC at a fixed 
transverse field T. We briefly mention the main ingre- 
dients used. We used 20 walkers, and performed im- 
portance sampling using, for each T, the optimized trial 
Boltzmann-like wavefunction ip^ opt " > discussed in Sec. Mil 
Branching was performed every 10 iterations for large 
values of T. For T < 0.25, however, weight instabilities 
are so severe that one needs to perform branching at ev- 
ery iteration. Finally, we made use of a continuous-time 
approach, sampling directly the probability of generat- 
ing an off-diagonal move with a Poisson's process^ In 
Fig. H we plot the GFMC results obtained for the 80 x 80 
random Ising model instance used in Ref. 0, for several 
fixed values of the transverse field T. The middle panel 
shows the GFMC estimate of the total energy (per site) 
for several values of T, compared to the variational re- 
sults, £m° ltz \ discussed in the previous section. The in- 
set allows to appreciate the differences between the two 
results in the small T region, which are invisible on the 
scale of the main plot. In the bottom panel we report 
several data regarding the residual diagonal energy. A 
small technical point is here in order. GFMC can calcu- 
late directly total energy estimates (or averages of other 
observables which commute with H) while averages of op- 
erators that do not commute with the Hamiltonian H are 
less straightforward to obtainiSS In particular, what one 
can simply evaluate is the so-called mixed averagei^^ of 
the potential energy H c i, 

\Hcl) w - ~. ~i i ~i v— > (!0j 

{Vt\Wgs) 

This is the estimator labeled "mixed" in Fig. The true 
expectation value we want is, instead, {iPgs\H c i\'4'gs}, 
which might be poorly approximated by the mixed av- 
erage if the trial wavefunction is poor. A simple partial 
cure to this drawback, is to include the so-called Ceperley 
correction^: {ipGs\H c i\ip G s) ~ 2(H cl ) w - (H c i) T where 
(H c i)t = {tpT\H c i\ipT) is the variational estimate ob- 
tained in the previous section. The results of the latter 
approximation for the residual energy are labeled "Ceper- 
ley" in Fig. 01 and are seen to be consistently lower than 



V. GFMC QUANTUM ANNEALING 

We finally present the results of a GFMC-based QA 
approach, where the transverse field T is decreased step- 
wise during the simulation, while, at the same time, 
the importance sampling Boltzmann-like wavefunction is 
changed according to the corresponding value of the vari- 
ational parameter /3 opt (T). 

As a benchmark, we will compare GFMC-QA out- 
comes with the Path-Integral Monte Carlo quantum an- 
nealing (PIMC-QA) and class ical simulated annealing 
(CA) results described in Refs. Ifil7t We reduce the cou- 
pling T in Eq. [T] at each Monte Carlo step (MCS) in 
a linear way: We start from an initial large enough 
value of the transverse field, To = 2.5, and then we set 
r„ = r (l - n/r) during the n-th MCS (0 < n < r). t is 
the total annealing-time (see Fig. ^) measured as the to- 
tal number of MCS performed by the algorithm. We used 
M = 20 walkers, and performed branching at each MCS, 
because the low-r region is affected by severe weight in- 
stabilities which would otherwise make the algorithm un- 
stable (for the initial, large-r, part of the annealing one 
could consider branching less often, as weights are well 
under control; this makes a negligible difference). How- 
ever, even with this very conservative choice, the weights 
go sometimes completely wild if a T of order 10~ 7 is 
reached (i.e. for long annealing-times r): As a conse- 
quence, we decided to cut-off the T n annealing schedule 
in such a way that the final T is w 10~ 4 and not 0. This 
finally guarantees a good weight stability. (Whenever it 
was possible to perform annealings with smaller cut-offs 
on r, and approximately the same slope 1/r, we checked 
that the results obtained are not very sensitive to the cut- 
off chosen.) For each value of r n the trial wavefunction 
used is the Boltzmann-like one, defined in Eq. 1151 with 
a variational parameter /3 pt(r n ), which corresponds to 
the instantaneous optimal value. (Practically, we used for 
/3(r) the fitting function shown in Fig. 01 upper panel.) 

Fig. El shows the best residual energy per spin ever 
reached during the annealing simulation, for several val- 
ues of t, averaged over 10 independent repetitions of 
the whole annealing process (due to computer time- 
limitations, a single run is shown for the largest, r > 10 8 , 
annealings). For comparison, the CA and PIMC-QA 
data obtained in Ref. kj are also shown. Notice first that 
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the r axes of the three calculations are completely un- 
related: The GFMC r is measured in units in which a 
MCS is just a single spin-flip, while MCS for the CA and 
PIMC-QA are intended as sweeps of the entire lattice of 
N spins (including all the P = 20 Trotter slices, for the 
PIMC case). For this reason, we also present the CA and 
PIMC-QA data in a shifted time axis where r is multi- 
plied by a factor N = 80 2 (rightmost curves). Although 
the GFMC-QA data are strictly below both the CA and 
the PIMC-QA data, on the same per-spin time- unit (i.e., 
compared to the shifted CA and PIMC-QA data), it is 
clear that the GFMC-slope is still worse than that of 
PIMC-QA, and indeed surprisingly similar to CA. More- 
over, the CPU-time needed for a single spin-flip in GFMC 
is much larger than the corresponding single-spin move 
in CA or PIMC-QA (each GFMC move costs of order N 
operations). 

Let us pause to consider the similarity between the 
CA and the GFMC-QA slopes that Fig. El suggests. This 
similarity must be somehow related to the fact that we 
have used, as importance wavefunction for the GFMC, 
a Boltzmann-like wavefunction, i/)t(x) oc e~^/ 2 ^ EcltyX \ 
At the level of a plain variational Monte Carlo (VMC), 
we already pointed out that such a choice of wavefunc- 



tion amounts to sampling \tpif (x)\ 2 — e P^W, and 



b^\x)\ 2 = e-Z^'W 
is therefore totally equivalent to a classical Metropolis 
Monte Carlo at temperature T = 1/(3. If, during the 
GFMC simulation, we neglect the weights associated to 
the walkers (as well as the associated branching process), 
we will be carrying over a completely classical simulation 
where the generated configurations are distributed (see 
Eq. ITT|) according to 



Gx'. 



^t\ x ') = Te -(fi opt (T)/2)[E cl (x')-E cl {x)\ 



As a consequence, the Markov process in x-space will 
obey a classical detailed balance condition 

Px',x e-^-M = px x , e -PoMx') (17) 

In other words, a GFMC-QA without weights would be 
just a computationally heavy way of doing a classical 
Simulated Annealing with a peculiar form of annealing 
of the effective temperature p op t(T) (notice, in passing, 
that such an optimal effective temperature never gets too 
low, since (3 op t saturates to around (3 op t ~ 2 for low T). 
Quantum Mechanics enters, therefore, only through the 
weights that the GFMC carries over (and the unavoidable 
branching process which makes the multiplicative process 
of weight updating numerically stable) . Evidently, such 
a weight updating is in the present disordered case not 
sufficiently strong and effective as to make the result- 
ing averages really different from the underlying classical 
physics governing the Markov chain in x-space, and the 
resulting GFMC-QA data are rather close to the CA ones 
(although they are much more computer-time demand- 
ing). 



VI. DISCUSSION AND CONCLUSIONS 

In this paper we have investigated the practical feasi- 
bility of a Green's function Monte Carlo (GFMC) as a 
tool for performing Quantum Annealing (QA). As a nat- 
ural test case, we have concentrated our attention on the 
two-dimensional Edwards- Anderson Ising model in trans- 
verse field, which was studied in Refs&i using PIMC-QA 
as well as standard CA. 

We identified the choice of the importance-sampling 
trial wavefunction, a necessary ingredient in any GFMC, 
as the crucial step - as well as the weak point - of a 
GFMC based QA (GFMC-QA). In particular, we found 
that the simplest mean-field wavefunction (analogous in 
many respects to the Weiss theory of ferromagnetism) is 
completely useless: Its optimization requires finding op- 
timal local fields hi (or magnetizations rrij = tanh(/ii) G 
(—1, 1)) describing the single-site wavefunctions, a com- 
plex problem with many minima which is as difficult as 
the original optimization problem, i.e., finding a clas- 
sical configuration of spins Si = ±1 optimizing the 
Edwards-Anderson classical term E c i = — JijSiSj. 
Using, instead, a simpler Boltzmann-like trial-function 
tPt{x) oc e~^/ 2 ) E ci( x ) (with f3 a tunable parameter which 
one can optimize variationally) , the resulting GFMC is 
feasible, but the corresponding residual energy results are 
disappointingly close - in magnitude and in slope, when 
considered as a function of the annealing time r - to 
those found by a standard classical simulated annealing 
(CA) (computationally much cheaper). We can rational- 
ize this finding with the inability of the GFMC algorithm, 
in the present disordered context, to properly implement 
the quantum mechanics inherent in the weights that the 
walkers carry with them. 

At the technical level, we also found that the choice of 
the number of walkers and the length interval between 
two successive branching events, can strongly affect the 
GFMC-QA performance, notably the stability of the al- 
gorithm. 

The crucial theoretical question is, therefore: How can 
one find good trial variational wavefunctions which de- 
scribe well enough the small-r glassy phase of an Ising 
spin-glass? This is, quite evidently, a highly non-trivial 
task. Taking inspiration from the existing literature on 
quantum models without disorder, one might think of in- 
troducing pair-correlations into the trial wavefunction - 
for instance, by means of spin-spin Jastrow factors, cither 
at nearest-neighbor or at longer range - as usually done 
in the framework of correlated lattice models^ and of 
electronic structure calculationsi22i2£ Unfortunately, for 
a quantum spin-glass, due to frustration and disorder, the 
form of such pair-correlations is far from obvious. More- 
over, whenever a large number of variational parameters 
in the trial-function is required, very advanced minimiza- 
tion techniques, such as those discussed in Ref. |23, are 
mandatory. This kind of computational schemes, how- 
ever, have been successfully tested only in equilibrium 
simulations of ordered systems, while our GFMC-QA 
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should cope, inherently with a non-equilibrium dynamics 
in a disordered system, a highly non-trivial step forward. 

We conclude that, at present, without a serious effort 
in constructing reliable importance sampling variational 
wavefunctions for a quantum glass, GFMC-QA is not a 
true competitor of PIMC-QA. 
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